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Using the nonequilibrium theory of superconductivity with the tunnel Hamiltonian, we consider 
a mesoscopic NISINISIN heterostructure, i.e., a structure consisting of five intermittent normal- 
metal (N) and superconducting (S) regions separated by insulating tunnel barriers (I). Applying 
the bias voltage between the outer normal electrodes one can drive the central N island very far 
from equilibrium. Depending on the resistance ratio of outer and inner tunnel junctions, one can 
realize either effective electron cooling in the central N island or create highly nonequilibrium energy 
distributions of electrons in both S and N islands. These distributions exhibit multiple peaks at 
a distance of integer multiples of the superconducting chemical potential. In the latter case the 
superconducting gap in the S islands is strongly suppressed as compared to its equilibrium value. 

PACS numbers: 73.23.-b, 74.78.-w, 74.45. +c 



I. INTRODUCTION 

Mesoscopic electronic applications typically rely on 
phenomena which show best when the electrons in small 
wires are cooled to very low temperatures, ideally to the 
range of 10 mK. In this regime the crystal lattice is very 
weakly coupled to the electron system, and electron cool- 
ing via the lattice becomes difficult. An alternative ap- 
proach is then to directly cool the electrons. This can 
be achieved by placing superconducting (S) contacts via 
insulating (I) barriers to the normal-metal (N) or super- 
conducting (S') wire whose electrons are to be cooled^^. 
By applying a voltage over such SINIS/SIS'IS coolers, it 
has been shown that one can cool electrons well below 100 
mK with these structures^, even when the lattice remains 
at a few hundred mK, or to enhance the superconductiv- 
ity in the middle S' island£*§i Optimally, such cooling 
should take the electron temperature to a few mK, a limit 
which is hardly reached in mesoscopic systems via other 
known means. 

With this type of nonequilibrium cooling, the concept 
of the electron temperature is not always well defined 4 , 
but one has to rather describe the full electron energy 
distribution function 8,9 . In this case cooling corresponds 
to the sharpening of this distribution function, essentially 
removing the high-energy excitations. 

One of the features limiting the performance of such SI- 
NIS coolers is the fact that the poor heat conductivity of 
the superconductors makes them inefficient reservoirs^. 
An additional pair of normal-metal electrodes attached to 
superconducting electrodes of a SINIS cooler can improve 
the relaxation and enhance the cooling characteristics of 
the device^. In this paper, we consider the effects of ex- 
tra N electrodes of this type attached to a generic SINIS 
structure. The superconducting pieces are now assumed 
small enough so that they can be driven out of equilib- 
rium by applying a bias voltage between the two normal 



electrodes. It is this arrangement of a NISIN'ISIN mul- 
tiple heterostructure which is the subject of the present 
study. Its novel feature as compared to the traditional 
SINIS structure with bulk S electrodes is that nonequilib- 
rium is now induced in all inner islands of the structure, 
which, in turn, strongly enhances nonequilibrium effects 
both in the central N and in the adjacent S islands. The 
resulting distributions in each island can be inspected by 
transverse probe tunnel junctions^. 

The microscopic nonequilibrium theory of double- 
barrier SINIS junctions is based on time-dependent 
Keldysh Green function formalism (see Ref. ^ and ref- 
erences therein). In the present paper we extend the 
theory such that nonequilibrium distributions in super- 
conducting islands are also allowed. The major modifi- 
cation is that the chemical potentials of the two super- 
conducting islands cannot be chosen zero simultaneously 
since these islands have different potentials determined 
by the bias voltage. This is equivalent to a time depen- 
dence of the order parameters imposed by different time- 
dependent order-parameter phases. We restrict ourselves 
to a tunneling Hamiltonian approximation which effec- 
tively makes the problem spatially independent within 
each superconducting or normal island. 



II. MODEL 

The system we study is shown in Fig. ^ In nonequi- 
librium, each of the islands has a separate energy distri- 
bution and, according to our model based on the tunnel 
Hamiltonian, they are independent of both coordinates 
and directions of the momenta. This implies that each 
region is in a diffusive regime when the momentum direc- 
tion dependence is averaged out within the first approx- 
imation in £/L where £ is the impurity mean free path 
and L is the length of the contacting region. In addi- 
tion, one has to assume that the intrinsic normal-state 
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FIG. 1: NISINISIN structure: The central normal island (N) 
is placed between two superconducting islands (S) which are 
connected to two outer N electrodes serving as external leads. 
Each connection is realized through a tunnel contact with a 
resistance Ri^. 



resistance r of each island is much smaller than any of 
the tunnel resistances R to satisfy the condition that the 
potential variation inside each region is smaller than its 
drop across the tunnel barrier. 

A characteristic rate for tunneling from region 2 into 
region 1 to be defined later is 7712 = [4^ie 2 f2ii?] , where 
v\ is the normal-state density of states (DOS) at the 
Fermi level in region 1 and Oi is the volume of conduc- 
tor 1. Competing with the inelastic relaxation, this rate 
determines whether the injection of new quasiparticles 
into the system is fast enough to drive the system out 
of equilibrium. The inelastic processes which uphold the 
thermal Fermi distribution can be neglected if 77 1/ Ti ne i 
which is equivalent to 



ve 2 QR/Ti r 



(14 /t 



2 ) 
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r<i 



(i) 



where T = (N m R/R Q )(A/A). Here A stands for the 
average cross-section area of the conductor, L is its length 
such that = AL, A is the area of the junction, £- me i = 
V-Dfinci, and D = vf£/3 is the diffusion coefficient. We 
also introduce the quantum resistance Rq = irh/2e 2 , and 
the number of conducting modes N m ~ p 2 F A/h 2 in the 
area of the junction. 

In our work, we assume that the electron energy dis- 
tribution within each island is homogeneous, which also 
implies a homogeneous potential and a constant effective 
temperature. Analysis of kinetic equations shows that 
this is achieved if r <C R or 



aA/L ~ uDe 2 A/L > 1/R 



which gives 



DL~ 2 > {ve 2 RViy 



V 



This puts a restriction L <C tT. The limit £ <C L holds 
if r 3> 1. Our model is thus applicable provided 1 <C 
L/£ <C r. Condition (JTJ of small inelastic interaction 
reads £? nel > TL£. 



III. FORMALISM 



A. Transport equation 



We use the standard Keldysh Green function formal- 
ism, see for example Ref. [j2 We denote the Nambu- 



space matrix operator 

d v 2 





A* 



G- 1 =t 3 t-- — -n + H, H = 
or 2m 

It acts on matrix Green functions 

/ • 1 G F 

( ' - 1 -Ft g 



where G stands for either retarded (advanced), G R ^ A \ 
or Keldysh, G K , Green function. For tunnelling between 
two superconductors (or normal metals) 1 and 2, the self 
energy in region 1 is St(1) = M7i2fl(2) provided the semi- 
classical Green functions 



9 



^G(p- 
ni 



k/2,p-k/2) 



integrated over the energy £ p = p 2 /2m — Ep, are inde- 
pendent of the directions of momentum p. The tunneling 
rate is defined as 



(2) 



Here T p (i is the tunnel matrix element and R12 is the 
tunnel resistance of the contact between regions 1 and 2. 
Note that the self-energy in region 2 is &r(2) = i"?2is(l) 

where 7721 = nv^ 1 (|Tp,q| 2 ) = [4z/ 2 e 2 fi 2 i2i2] ~ ■ If re- 
gion 1 has tunnel contacts to two other regions 2 and 3, 
the self energy is a sum 

E T (1) = ir)\ig{2) + w?i 3 <7( 3 ) 

where 7712 = [4^ 1 e 2 r2 1 i? 12 ] 1 and i] 13 = [4^ 1 e 2 r2 1 i? 13 ] . 

Including the tunnelling self-energy into the total self- 
energy that contains both elastic and inelastic processes 
we can write equations for the retarded, advanced, and 
Keldysh Green functions: 

(g- 1 - £*< A >) o g r w = is( Xl - X2 ) , 

G- 1 - ± R ) o G K — Ti K o G A = . 



The product t K o G A is a convolution over frequencies 
and momenta of the type 



t K o 



G = 



ei,e e,e 2 ^ 



Applying the inverse operator to G R W and G K from 
the right and subtracting the obtained equations from 
the equations above we find the transport-like equations 



G 



K 



,£2 1 "61,63^2^*3 

, G K_ G K 



H o G K — G K o H = X* , (3) 



v^kG^-erfsGf^+Gf^e 

+H o G R W - G R W o H = l R [fJ , (4) 
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where the collision integrals in each region are 
i* ue2 = t R oG K -G K ot A -G R ot K + t K oG A 

jR(A) = £.R(A) Q qR(A) _ qR(A) q £R{A) _ 



We also define 



h h 

4 



-4 h 



Finally, the tunnel collision integral for region 1 in con- 
tact with region 2 takes the form 

#(1) = % 2 [g R {2)o g K (l) -g R (l)o g K {2) 

+g K (2)og A (l)-g K (l)og A {2)] . (5) 

If the distribution function is time independent, the 
Keldysh Green function has the form 

9e u e 2 = 9? u e 2 (A,e 2 +^3/2,62) ~ (A,6i + hf2,ei) <jt u e 2 > 

(6) 

where the odd- and even-in-e distribution functions /i 
and ji are defined such that 



A 



-n(e) + ra(-e) , f 2 , e = 1 - n(-e) - n(e) 



The actual quasiparticle energy distribution function n 
satisfies fi <£ + f 2 , e = 1 — 2n(e). 

For spatially uniform distributions the kinetic equa- 
tions in each conducting region is obtained after averag- 
ing Eqs. Js3, over the momentum directions. The elas- 
tic collision integrals average out, and the kinetic equa- 
tions for diagonal in e components become 



Tr 



Tr 



fK 



I" - I* 



A 



Tr (f 3 M £ , e ) = Tr (f 3 [l K ~ (/£ - I A ) A, 



(7) 
(8) 



which is the conservation of charge Q(l) = eiV(l)f2i in 
the volume Qi of region 1: 



gg(l) 
0t 



= 1(1) 



(9) 



Here -/V(l) is the particle density in region 1. The current 
that flows into region 1 is 



, de dQr, 
I(l) = -ieu Fl Qi I ^^Tr 



r 3 /f(l) 



(10) 



The inelastic collision integral drops out because it con- 
serves the particle number. 

Multiplying Eq. by e and taking the trace of it one 
can obtain the balance of the energy £ = flE in the form 



at 

Here the energy density is 



de d 3 p |A| 2 



and the energy current into region 1 is 



Mi) 



de dfl 



4 47T 



p Tr 



(e- eV >ifa)I A (l) 



(11) 

The collision integral in Eq. contains both tunnel 



and electron-phonon contributions, / A '(l) 



- P h 



(1) 



(1). The electron-electron interactions drop out be- 
cause of the energy conservation. The energy flow into 
region 1 can be separated into two parts. One part con- 
taining l£L ph is the energy exchange with the heat bath 
(phonons). The other part contains the tunnel contribu- 
tion 1^ (1) and is the energy current into region 1 through 
the tunnel contact. 



The collision integrals here include only tunnel and in- 
elastic contributions I — It + And , where inelastic pro- 
cesses include electron-phonon and electron-electron in- 
teractions, /inei = ie-ph + ^e— e- The matrix M is 



M = [Hog K -g K oH 



H o g R - g R oH)f l + h (Hog A -g A o h) 



B. Charge and energy tunnel currents 

After taking the trace and integrating over energy and 
momentum, Eq. © in region 1 yields 



a_ 

at 



-IV\ 



de f d^ dn p ^ £ K 
2 I m An 



de df2 F 



-Tr 



IV. KINETIC EQUATIONS IN HYBRID 
STRUCTURES 

A. Nonzero superconducting chemical potential 

If a hybrid structure containing more than one super- 
conducting island is voltage biased, at least one of the 
superconductors will have a nonzero chemical potential 
with a time dependent phase of the order parameter. If 
the external voltage is constant and the resulting state is 
stationary the phase is linear in time, \ = while 
the magnitude of the order parameter is time indepen- 
dent. Consider two sets of regular Green functions. One 
set of functions taken in the basis where the 

order parameter phase varies in time. The other set be- 
longs to the same order parameter magnitude but has a 
phase constant in time (zero for simplicity). However, 
the energies of the Green- functions g R ^ A \ g R ( A \ f R ( A \ 
and /t-R(-A) are shifted by different amounts proportional 
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to the chemical potential in the previous representation. 
The relations between these two sets, the Green func- 
tions gf\f^ and f^i and the steady-state functions with 
shifted energies, can be established from the Eilenberger 
equations. Since A = |A|e +l2AIS ', i.e., 

A w - \A\2tt5(lu + 2fi s ) , A* = \A\2ttS(lu - 2 M5 ) , 

one observes that Eqs. 10} are satisfied by 

and 



(13) 



R(A) 



The functions g. 
isfy the steady-state Eilenberger equation 



-gf^ and 



-2ef^ + 2\A\g^=I 2>e 
supplemented with the normalization 



R(A) 



f! R(A) sat- 



(14) 



-MA) 



jR(A) 



(15) 



Therefore, we find from Eq. JBJ) for a nonzero chemical 
potential 



g* 6 , = - (/ M + / a , e ) 2^ - e') , (16) 
- (A,e - /a, e ) 2^(e - . (17) 



The off-diagonal Keldysh Green function in Eq. JBJ be- 
comes 



f A 



(/i.e + /2, e )/£u s ] 2tt5(o; + 2 Ms ) . (18) 



B. Tunnel collision integrals 



Assume that region 1 is a superconductor and region 
2 is a normal metal. With Eqs. {TSJ), ifTflJl. and lfT7)l. the 
matrix elements of the tunnel collision integral Eq. JSJ in 
the superconducting region S give 



Tr 



if (5) - 4t VsN {[ge+^iS) + g^ s (S)} \f u (S) - f he (N)] + [g e+fls (S) ~ fl e - M (S)] [/ 2 , e (5) - / 2 , e (iV)]} , (19) 



Tr 



f 3 7#(S) = 4i7/ SJ v {[ 5e+ „ s (S) " fc-wOS)] - hAN)} + [ge +t , s (S) + ge-^S)] [/ 2 , e (5) - / 2 , e (JV)]} . (20) 



Components of the retarded and advanced tunnel collision integrals are 

I« {A \S) = ±i VSN f R W(S) , -I™ A \S) = ±i VSN f^ A \S) , I^ A) (S) = I* A \S) = . 
Both Keldysh and R(A) tunnel collision integrals in the normal region N are coupled to those in the region S by 

vsnI(n) = -msi(S). 

Using Eqs. |p5J and JTEJl, we find 

Tr (M e , e ) = 2|A|F e+Als (/ Xj£ - /i, e+2w + / 2)E + /a.e+a^) - 2|A|F e * MS (/i, e _ 2/lfl - h,e + /a,e + /a,e-a Ms ) , (21) 
Tr (r 3 Me,e) = 2|A|F e+Als (/ 1)£ - /i, e +2 Ms + / 2)E + /a,«+a M ) + 2|A|^ IS (/i, e _a Ms - / 1)£ + /a,« + /a, e -a M ) • ( 22 ) 



r 



After inserting Eqs. (dJ), $2® and (EH into Eqs. © 
and JHJ they yield the final kinetic equations to be solved 
for the distribution functions. 

In the equations above, we introduce g e = 
(g?-g*)/2 = Reg? and also F e = (/* + f£) /2 = 
iIm/ £ . The functions g e and i 7 ^ are even in e. With the 
account of tunnel and electron-phonon interactions, they 
can be found from the steady-state Eilenberger equations 
(|14fl . i|15|) . These equations can be more easily solved 
when the inelastic interaction is small. If region 2 is nor- 



mal, one finds 



g, 



R(A) 



fR(A) 

where ~f = t]sn- 



± 



= ± 



e ± 17 



(e±^7)-|A| 2 
|A| 

(e±z 7 ) 2 -|Ap 



(23) 
(24) 
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For a N1IS1INIS2IN2 structure, each conductor Si, N, 
or S2 has tunnel contacts with two other conductors, thus 
the tunnel integrals are sums of the contributions from 
two regions. For example, if the inelastic interaction is 
small, the depairing rate in Eqs. I|23l) and (|24|1 is j(Si) = 

VStN + VS.N,- 



Current conservation 



For a NS junction, the tunnel current Eq. (|10() into the 
superconductor S from a normal region N or Ni,2 is 



I(N -> S) 



1 



deg e (S) [/i, e - MS (S) - fi,e+„ s (S) - A,e-» S (N) + fi, e+lia (N) 



— 00 



Electric currents through both junctions are balanced when I(N — > Si) + I(Ni — > Si) = 0. 

The energy current flowing into the superconducting island through each junction has the form 



1 



deeg^S) ([/i, e _ MS (5) + /i, e+MS (5) - fi, c -^(N) - f 



1,6+/*, 



f (iV)] 



(25) 



+ [/2 )e - M (5) - / 2 , e +i. s (5) - f2, e - Ma (N) + h,e+»s(N)]) ~ H N - ^)[Ms/e + p s ] • 

I 



(26) 



The energy conservation follows from the kinetic equa- 
tions and is therefore an abundant condition. However, 
in the quasi-equilibrium limit, when the electron-electron 
interaction in each island dominates over the tunnel injec- 
tion, the distribution functions are nearly thermal with 
certain temperatures specific for each region. In this case, 
the kinetic equations do not need to be solved explicitly, 
but the temperatures are found by requiring the conser- 
vation of the energy current. 

The energy current Is(S — * N) flowing into the cen- 
tral normal island through each junction is obtained 
from Eq. i|26|) by changing the sign and replacing ips 
with ipN so that the difference between the energy cur- 
rents from N to S calculated in two opposite directions 
is I e (N S) + I e (S -> N) = I{N -► S)V NS , where 
Vns — <Pn ~ fs is the voltage across the junction. This 
difference of energy currents is compensated by the work 
produced by the voltage source and does not lead to the 
change in the energy of the superconductor. When the 
energy currents into a superconducting island are bal- 
anced, the terms with [is/e+<Ps cancel due to the electric 
current conservation. In the normal island, however, only 
the terms with ip^ drop out since the chemical potentials 
^5 on both sides of the island are different. 



D. Self-consistency equation and charge neutrality 

The magnitude of the order parameter is found from 
Eq. iJTSJl. It is a real quantity, therefore, 



A 



Re jf) [fx 



e+MS ~r 71,6— (is 



de 

T 



(27) 



where E c is the BCS cut-off energy. The self-consistency 
of the order parameter requires vanishing of its imaginary 
part: 



/oo 
Fe(fl, 
-OO 



fx, 



which results in 

TrM de 



Tr 



2,e+/i S 



h,€-ns) de = , 
(28) 



de = 



(29) 



Equation (|28|l or the second condition in Eq. Ij29() to- 
gether with the kinetic equation Eq. JSJ is equivalent to 

jTr(r$i K j de = 0. Since the integrals over the inelastic 
electron-phonon and electron-electron collision integrals 
vanish, this is the condition of current conservation: 



Tr ( f 3 /|F ) de = 



(30) 



which implies through Eq. (|10() that the sum of currents 
flowing into an island is zero. 

Charge density can be obtained by calculating Trg K 
from Eqs. JTSJ, Q17[l. Making shift of e in the integral we 
obtain 



N = No - 2v (eip + us) 

9e (/l,e+/ig — fl,e— fj,s 



f: 



2,6+Ms 



de 



2,6-Ais^ 



Here we take into account the divergence of the integral 
for large e. Charge neutrality requires that N = N . The 
electric potential (fs of a superconductor is thus coupled 
to the chemical potential through 



e$ = 



9e (fi, 



e+MS 



/l,6 



2,e+Ms 



f2,e—ns) ^ 7 

(31) 
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FIG. 2: (a) Chemical and electrostatic potentials of the su- 
perconducting island as a function of bias voltage, (b) I — V 
curves for different p. For clarity, in the case p = 1000 the 
potentials have been multiplied by the factor 100 and the cur- 
rents by the factor 1000. 

where e$ = ps + eips- 

To summarize, for a system of n superconducting and 
m = n + 1 normal pieces we have n + m pairs of the 
distribution functions /i and f 2 and n parameters pis- 
For these unknowns, we have 2(n + m) functional kinetic 
equations Eqs. J7J), JSJ) and n conditions Eq. (|30|l of cur- 
rent conservation for each superconducting island. Note 
that the current balance of the type of Eq. (|30fl in each 
normal island is satisfied automatically due to the kinetic 
equations that have M = in any normal conductor. 

V. RESULTS 

To further simplify our model, we limit the discussion 
to a case symmetric with respect to the central normal is- 
land, i.e. Rn x Si = Rn 2 S 2 = Rout, RnS ± — RnS 2 = Rin, 
| Ai| = |A 2 | = A, and cp Nl = -<fN 2 , MSi = ~^s 2 = Ms- 
Therefore /i(JVi) = fi(N 2 ) and f 2 (N 1 ) = -f 2 (N 2 ). 
Moreover, when applied to kinetic equations in the nor- 
mal region N, the symmetry yields A (Si) = —f 2 (S 2 ) = 
MS) , A (Si) = /i(S 2 ) = /i(S) and f 2 (N) = 0. 

A. Nonequilibrium state 

1. Distribution functions 

Under the conditions formulated in Sect. II when the 
inelastic relaxation can be neglected, r\ T^_ e , the 
kinetic equations J7J and © were solved numerically 
for nonequilibrium distributions in both normal and su- 
perconducting islands. Throughout numerical computa- 
tions we use the value of the pair-breaking parameter 
A/7 = 1000. For nonequilibrium states, the bath tem- 
perature was fixed at a value much lower than the gap, 
Tbath = (0.1/1.764) A. This corresponds to T = 0.1 T c if 
the gap coincides with the zero-temperature BCS value 
A = 1.764 T c . 

The distribution functions in both N and S islands are 
in turn determined by the chemical potential of supercon- 
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FIG. 3: Nonequilibrium normal-island distributions for p = 
1000. We put v = eV/A here and in Figs. EHZI below. 



ducting islands /is which was calculated self-consistently 
using the current conservation as discussed in the previ- 
ous section. Chemical potential and the current -voltage 
curves are shown in Fig. [5| 

We find three qualitatively distinct cases characterized 
by the ratio p = R ou t/Rin- For the ratio as large as 
p = 1000, the distribution in the normal island, njy, is 
shown in Fig.|21for several values of V . Under biasing V, 
the distribution is driven into nonequilibrium. At first, 
this is seen only as a slight deviation from the Fermi 
function np but for voltages V/A > 1, the distribution 
njv assumes the characteristic step-like profile. This be- 
havior of npf can be explained as follows. For very high 
ratios p, the superconducting chemical potential is small: 
Ps/eV <C 1 [compare with Fig. 0(a)]. The kinetic equa- 
tions 0, © yield then f^N) « A(S) w AW). Since 
/ 2 (A) = 0, one has 2n N = 1 - fi(N) = 1 - /i(Aq). 
The external leads Ai and A2 are in equilibrium such 
that /i(Ai) = -n F (e - eV) + n F (-e - eV). Thus, 
tin = [«f(c — eV) + n F {e + eV)]/2. Since f 2 (S) is small 
due to the high resistivity ratio, one has also rig ~ 
for the superconducting island. 

The distributions for the ratios p — 1 and p = 1/1000 
are shown in Figs.E]and[S] They also change their shape 
strongly above a certain voltage. The changes correspond 
to the sharp rise in the electric current through the junc- 
tion, Fig. 0b). The origin of this behavior is discussed 
in connection with the charge imbalance, see next sub- 
section. The distribution functions of the N island show 
a cooling behavior: The distribution becomes very steep, 
thus corresponding to a low effective electron temper- 
ature at such bias voltages when the chemical potential 
difference across the superconductor/normal-island junc- 
tion is \p,s\ ~ A. One can define an effective temperature 
through^ 



k B T cS = ^^J^ [l-f u {N)]ede . (32) 

keeping in mind that f 2 (N) = for the center N island. 
This definition gives the actual electron temperature in 
(quasi)equilibrium. The effective temperatures have min- 
ima for p < 1 as seen in Fig. 
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FIG. 4: Nonequilibrium normal metal distributions for p = 1. 




FIG. 6: Effective temperature T e fi of Eq. 1521 . 
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FIG. 5: Nonequilibrium normal metal distributions for p = 
1/1000. 




FIG. 7: Nonequilibrium superconductor distributions for var- 
ious bias voltages when p = 1. 



For a good contact between the outer normal elec- 
trodes and the superconductors S, i.e. when p = 1/1000, 
there is no deviation between the superconductor distri- 
bution function ns and that in the normal reservoir, a 
Fermi function with ip^ — V/2. This suggests that, for 
small p, the considered structure is similar to a SINIS sys- 
tem with superconducting reservoirs. The distribution in 
the central N island for low p is shown in Fig. [SJ It re- 
sembles the distribution found for a SINIS structure^. 

When the ratio p increases, the distribution in the su- 
perconductors, ns, deviates from the Fermi function as 
shown in Fig. \7\ Nonequilibrium distribution in the su- 
perconducting regions on both sides of the central normal 
island drives the state in the N island yet further from 
equilibrium, see Fig. 0] For larger p the cooling behav- 
ior becomes less pronounced and finally disappears, see 
Fig. 

For low and intermediate values of p, we observe novel 
features of highly nonequilibrium distributions both in 
the central normal and in the side superconducting is- 
lands. For small p — 1/1000, peaks in the energy distri- 
bution appear at energies ±e = A + ps (see Fig- - In 
addition to these, new peaks appear in Fig. 0] for larger 
p = 1 at energies ±e = — A+p,s, for voltages considerably 
exceeding A/e. Both sets of peaks come as a result of 
recursion from singularities in g t and F e at e = ±A in the 
kinetic equations. The new peaks at ±e = —A + ps are 
present as long as the distribution functions of the two 
reservoirs differ at the corresponding energies and appear 
as a result of suppression of the distribution function due 
to a large factor (A/ r y)F e ± fJ _ s in the sub-gap region in 



Eqs. (|21|1 and (|22|l . Thus, the requirement of the new 
peaks is roughly A — ps < eV/2 which can be fulfilled if 
ps 7^ — eV/2 as is the case for intermediate values of p. 



2. Charge imbalance 

As mentioned above, the distribution function of the 
center N island suffers a drastic change above a certain 
voltage. This change coincides with the upturn of the 
current as a function of the applied voltage in Fig. |2Ib) 
and is accompanied by a deviation of the chemical po- 
tential ps from the electric potential —etps in the adja- 
cent superconductor as determined by Eq. 131|) . In equi- 
librium, their difference $ = 0. In nonequilibrium, a 
difference between ps and —e(fs appears according to 
Fig. [2Ja). The singularities appear when the chemical 
potential difference between the superconductor and one 
of the contacting normal conductors approaches A. This 
corresponds to eV ~ 2A for a large mismatch between 
i?in and i?out or to eV/2 ~ 2 A for i?; n = i? ou t- To mea- 
sure potentials tps, a capacitive connection would be re- 
quired, in addition to the usual resistive connection only 
capable of detecting /15. 

For a good contact between the outer normal elec- 
trodes and the superconductors S, i.e. when p = 1/1000, 
there is no deviation between —eips and ps- This can be 
understood by considering the function ns, which coin- 
cides essentially with a Fermi function shifted by p>s ~ 
—eipN 1 = —eV/2. For a Fermi function, n e = 1 — n_ e , 
and the term in the brackets in Eq. I|31|l vanishes. 
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eV/A eV/A ( 



values with the dimensions of energy (like voltage, exci- 
tation energy, temperature) to the real gap magnitude A, 
we need to rescale the real voltage to its relative magni- 
tude. Conversion between the relative, V/A, and the real 
voltage normalized to the BCS gap, V/Aq, is provided by 
Fig. IHl^b) where the graphs are given for p = 1000 and 
p = 1. As seen from Fig. Etb), for p = 1 high relative 
voltages eV I A can be achieved for comparatively low ab- 
solute voltage values eVj A . 



FIG. 8: (a) Gap A as a function of bias voltage for two re- 
sistance ratios p = 1000 (dashed line) and p — 1 (solid line). 
Both quantities are normalized to zero-temperature BCS gap 
Ao. (b) Conversion between eV/ A and eV/ Ao for p = 1000 
(dashed line) and p = 1 (solid line) . 



3. Gap instability 

In a noncquilibrium state, the gap function A has to 
be calculated self-consistently using Eq. (|T7Jl . Employing 
the equation for the critical temperature T c , 



tanh ■ 



de 



2T r 



(33) 



one can exclude the interaction constant in favor of T c . 
For a low resistance ratio p = 1/1000, the gap does not 
change considerably, A ss A for all V, A = 1.764 T c be- 
ing the zero-temperature BCS gap. However, for a higher 
resistance ratio, the noncquilibrium energy gap is mod- 
ified dramatically as shown in Fig. |SJa). One observes 
a drastic change in the gap for voltages coinciding with 
those where the change in the distribution is seen. The 
energy gap becomes a multi-valued function which im- 
plies hysteretic behavior accompanied by jumps of A at 
the corresponding voltages. For a very poor contact be- 
tween the superconducting islands and the outer normal 
reservoirs, i.e. for high p when deviation from equilib- 
rium is the largest, the gap function jumps down to very 
small values and superconductivity is nearly destroyed. 
For a lower tunnel resistance ratio p, the gap decrease 
is not so huge and superconductivity is less suppressed. 
Note that with respect to the order parameter magni- 
tude, the bath temperature chosen for our calculations 
can be considered as zero. Indeed, the temperature was 
set much lower than A while the relevant energy scales 
for the distribution function are determined by the ap- 
plied voltage and by A itself. Thus the thermal effects 
on the gap are negligible. 

The predicted suppression of superconductivity in a 
nonequilibrium superconductor placed into a tunnel con- 
tact with a nonequilibrium normal-metal electrode con- 
trasts to the superconductivity enhancement observed 
in tunnel SIS'IS structures£*&£ where the nonequilibrium 
superconductor is in contact with equilibrium supercon- 
ducting electrodes. 

Note that since in our calculations we normalize all the 



4- Visualization 

The peaks in the distribution of N and S islands can 
be monitored by measuring the differential conductance 
of a probe tunnel SINIS or SIS'IS junction attached to 
the island in question. Let us consider a SINIS probe 
junction attached to the central N island as in Ref. 4. 
The distribution function in the N island is not modified 
by the measuring current if the tunnel resistance of the 
probe junction satisfies Rs p n 3> -Rin + -Rout- The current 
through the probe junction is 



I(N -» S P ) 



1 



9e(Sl 



eRspN j-c 
x [np(e) - n N (e + eVp/2)] de, 



where the bulk superconducting probe electrode Sp is as- 
sumed to be in equilibrium with a potential so that 
ps p = —etpp = —eVp/2 where Vp is the voltage be- 
tween the two probe electrodes. The energy gap Ap 
in the probe electrodes is assumed to have the magni- 
tude corresponding to the BCS value for T = Tbath- As 
?bath ~ 0.1 T c , the gap Ap is very close to A . In addi- 
tion, we set 7p = 7 for the depairing rate in the probe 
electrodes. The differential conductance becomes 



dl 



1 



2R Sl 



N 



dn N (e + eV P 2) 

9e [bp) 3 de . 34) 

de 



Due to the peaks in g e (Sp) at e = ±Ap, the differen- 
tial conductance should reproduce the peaks in the dis- 
tribution function at the probe voltages satisfying 
e±A P = Vp/2. 

The differential conductance dI(N — ► Sp)/dVp to- 
gether with the corresponding distributions in the central 
normal island for eV/A = 10 are shown in Fig. [§]for the 
ratio p = 1. The peaks in the distribution in Fig.[!Jb) are 
located at e = ±A - p s , he., e = 2.89 A and e = 4.89 A. 
Two more are located at e = ±A — 3ps or e — 10.7 A and 
e = 12.7 A (out of scale in Fig[3Jb)). Comparing these 
to the locations of the larger peaks in Fig.^a) and using 
the A-vs-Ao conversion curve of Fig. [Sf a) we see that 
the peaks in the distribution are indeed reproduced at 
the probe voltage eV P /A = 2e/A±2A P /A. The smaller 
peaks in Fig. Efa) refer to the less pronounced structure 
in the distribution not well-resolved in Fig. Elb). 
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(a) (b) (a) (b) 




FIG. 9: Differential conductance for the probe junction (a) 
and the corresponding quasiparticle energy distribution in the 
normal-metal island (b) when p = 1, eV/ A — 10. 



FIG. 10: Temperature in the normal-metal island as a func- 
tion of bias voltage for bath temperatures 0.05 A (a) and 0.2 A 
(b) when the system is quasiequilibrium. 



B. Quasi-equilibrium 

As discussed in Section II, the condition for full 
nonequilibrium is r\ r^, where rj is defined accord- 
ing to Eq. J3J). The inelastic relaxation may be fur- 
ther separated to relaxation caused by electron-electron 
and electron-phonon interactions with collision rates T~_} e 
and Tgl^, respectively. The experimental situation in 
nanoscale heterostructures^ corresponds often to the case 
t~_} c > t~ b . In particular, the limit of low coupling 
to the heat bath in S and/or N islands is frequently 
realized when the tunnel injection rate is intermediate 
between the electron-phonon and electron-electron re- 
laxation rates, t~ h <C rj <C 7^1 e . We refer to this 
case as quasi-equilibrium. While the near absence of 
electron-phonon interactions prohibits the quasiparticles 
from coupling to the lattice, the rate of electron-electron 
scattering is high enough for the quasiparticles to assume 
a Fermi distribution with certain electron temperature. 
We have studied the cooling performance of our NISIN- 
ISIN hcterostructure in the quasi-equilibrium limit look- 
ing at the electron temperature of the central N island. 

We note that the simple expressions for the regular 
Green functions of the form of Eqs. (|23|) and (|24[1 are 
not applicable in the strict sense when the inelastic re- 
laxation dominates. To find the exact expressions for the 
regular Green functions one has to solve the Eilenberger 
equations (|14|> , l|15|l with the proper inelastic collision in- 
tegrals. However, to simplify our problem, we model the 
pair breaking effects of inelastic relaxation by an effec- 
tive pair-breaking rate 7 in Eqs. (E3 and O in the same 
way as for the tunnel limit described in the previous sec- 
tions. This approximation is frequently used in practical 
calculations. Here we put A/7 = 1000 as above. 

Applying the tunnelling model shows that depending 
on the configuration of the quasiparticle traps, i.e. the 
ratio of outer and inner junction resistances p, effective 
cooling of the normal-metal island can be achieved as 
demonstrated also in Ref. 0. The temperatures of the 
central N island and of the contacting S island are shown 
in Figs. and ^2 respectively. The temperature of N 
island indeed displays a minimum below bath temper- 
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FIG. 11: Temperature in the superconductor as a function of 
bias voltage for bath temperatures 0.05 A (a) and 0.2 A (b) 
when the system is quasiequilibrium. 



ature when the ratio p is small. However, for large p, 
the temperature monotonously rises above the bath tem- 
perature with an increasing voltage. As can be deduced 
from Figs. ITUI and ITT1 the cooling effect is not attributed 
to the presence of two additional normal metal reser- 
voirs. On the contrary, smallest ratio p, which corre- 
sponds to the strongest cooling, is seen to lead to almost 
constant T$ — Tbath as it would be the case for pure su- 
perconducting reservoirs. Another remark concerns the 
cooling efficiency of a NISINISIN configuration as Tbath 
is lowered. Indeed, for Tbath = 0.2 A the temperature 
minimum for the depairing parameter 7 = A/1000 used 
for our calculations is roughly 0.063 A while the mini- 
mum is T min = 0.024 A for T ba th = 0.15A, and it is 
T min ~ 0.013 A for Tbath < 0.1 A. From the numerical 
results in the case 1 » p > 7/ A, we find that the min- 
imum achieved temperature follows T m i n /T c w 0.24p < = 
where £ « 0.5. For smaller i? ou t, the minimum tem- 
perature T m i n is determined by the inverse proximity ef- 
fect described by the depairing parameter 7 (see Ref. 0) . 
Combining both the superconductor heating due to a fi- 
nite i?out and the effects of depairing, we can write an 
approximate formula for the minimum temperature for 
relatively low bath temperatures, 



T min /T c = 2.5 ( 7 /A) 2/3 + 0.24P 1 / 2 



(35) 



For p <C f, the depairing rate 7 is limited by i? ou t, i.e. 



10 




0.5 1 1.5 2 2.5 3 3.5 4 

eV/A 




0.5 1 1.5 2 2.5 3 3.5 4 

eV/A 



3 
2.5 

2 
1.5 

1 

0.5 




(a) 



p=10" 5 




---p=10" 4 




p=10 J 




p=10" 2 




p=10" 1 




J 


i 



0.5 1 1.5 2 2.5 3 3.5 4 

eV/A 




0.5 1 1.5 2 2.5 3 3.5 4 

eV/A 



FIG. 12: Chemical potential of the superconductor as a func- 
tion of bias voltage for bath temperatures 0.05 A (a) and 0.2 A 
(b) when the system is in quasiequilibrium. 



FIG. 14: Differential conductance as a function of bias volt- 
age for bath temperatures 0.05 A (a) and 0.2 A (b) when the 
system is in quasiequilibrium. 
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FIG. 13: Electric current between S and N islands as a func- 
tion of bias voltage for bath temperatures 0.05 A (a) and 0.2 A 
(b) when the system is in quasiequilibrium. 



7 = l/(4v s e 2 £l S Rout) = (r s /2R out )E T , where E T = 
D/L 2 is the Thouless energy in the superconducting is- 
land with length L and r$ is its resistance. Substituting 
this into Eq. 135(1 , we find that the minimum temperature 
is optimized with 



pn6A[(E T /A)(r s /Rin)} 



4/7 



(36) 



Note, however, that Eq. l|3T)|l is valid provided rs/R- m <C 
p< 1. 

The sharp rise in T/v and T$ occurs generally around 
V = 2 A but for larger values of p, when also Tb a th — > 
T m j n , the upturn shifts towards higher voltages V. This 
is because the upturn is determined by the condition 
—ps ~ A rather than V/2 — A as can be seen by com- 
paring Figs. and ^2 to Fig. ^] For increasing bath 
temperature, though, this trend is smeared and disap- 
pears. The voltages corresponding to the temperature 
rise are also seen in the IV-curves in Fig. 1131 and in the 
differential conductance, Fig. ^] 



VI. DISCUSSION 

A. Electron cooling 

According to our results, the electron cooling is the 
most effective when the outer resistance is low p « 1. 



In fact, both the effective temperature and the distribu- 
tion function in the superconducting region almost coin- 
cide with those in the bath for such voltages that yield 
\p s \ < A especially for very low ratios of p. When the 
ratio p is low, a good contact between the inner super- 
conducting region and the outer electrode makes the dis- 
tributions in these two regions not so much different from 
each other, thus decreasing the role of the extra junction. 
This conclusion is valid only within the tunneling approx- 
imation. When the contact between the outer electrodes 
and superconducting islands are more transparent, cool- 
ing properties of the device are affected by the inverse 
proximity effect from the external normal leads. 

For larger ratios p, the extra junction prevents the 
state of the superconducting region from reaching equi- 
librium, thus reducing the cooling power of the entire 
structure. Moreover, this limit has another disadvan- 
tage as far as the cooling performance is concerned: For 
larger voltages when eV approaches 2A, one expects a 
suppression of superconductivity in the S regions down 
to lower values of A and thus the cooler would become 
even less effective. This suggests that the cooling per- 
formance of an NISINISIN structure cannot be improved 
essentially by an extra tunnel junction as compared to 
that of a simple SINIS structure. However, the presence 
of the quasiparticle traps helps to practically realize the 
superconducting reservoirs by thermalizing them quickly 
to an object with a high thermal conductance. 



B. Nonequilibrium distribution 

Nonequilibrium distribution formed in the supercon- 
ducting islands for high values of p results in yet stronger 
deviation from equilibrium in the central normal island. 
As seen from Figs. 0] and [S] the distribution function 
in the N island is characterized by peaks at energies 
±e = A ± p s , and also ±e = A ± 3/i s , etc., for voltages 
considerably exceeding A/e. These peaks are clearly vis- 
ible in the differential conductance of a probe tunnel SI- 
NIS junction made at the central normal island, Fig. [5] 
Simultaneously, nonequilibrium states in the supercon- 
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ducting region follow the gap which is strongly reduced 
as compared to its equilibrium BCS value A . The tran- 
sition into a nonequilibrium state is accompanied by a 
jump in the gap magnitude which leads to the jump in the 
relative voltage: high values of eV/A can be reached al- 
ready for not very large absolute values of voltage eV/Ao. 
This makes observation of the nonequilibrium states in 
N and S regions more easily accessible in experiments. 
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